A fully coupled system of generalized thermoelastic theory for semiconductor medium

This study presents a new mathematical framework for analyzing the behavior of semiconductor elastic materials subjected to an external magnetic field. The framework encompasses the interaction between plasma, thermal, and elastic waves. A novel, fully coupled mathematical model that describes the plasma thermoelastic behavior of semiconductor materials is derived. Our new model is applied to obtain the solution to Danilovskaya’s problem, which is formed from an isotropic homogeneous semiconductor material. The Laplace transform is utilized to get the solution in the frequency domain using a direct approach. Numerical methods are employed to calculate the inverse Laplace transform, enabling the determination of the solution in the physical domain. Graphical representations are utilized to depict the numerical outcomes of many physical fields, including temperature, stress, displacement, chemical potential, carrier density, and current carrier distributions. These representations are generated for different values of time and depth of the semiconductor material. Ultimately, we receive a comparison between our model and several earlier fundamental models, which is then graphically represented.


List of symbols α n
Discrepancy in deformation potential between the conduction and valence bands

Displacement components
In contemporary semiconductor technology, the shrinking size of devices poses challenges for thermal control.Consequently, it is imperative that we acquire a comprehensive comprehension of interconnected thermoelastic phenomena.The production of high-density integrated circuits (ICs) serves as a crucial illustration of this phenomenon.During the production process, specifically with methods like rapid thermal annealing (RTA), semiconductor wafers undergo sudden temperature fluctuations to activate dopants or rectify implantation damage.This process of thermal cycling results in significant temperature variations within the wafer.This results in substantial thermoelastic strains due to the mismatch in thermal expansion coefficients between the various layers of material.Applying these stresses can hurt the wafer in many ways, such as by warping it, making tiny cracks in it, and even causing major failures like delamination where the semiconductor meets other materials, like dielectrics or metal interconnects 1-3 .Plasma, thermal, and elastic models are becoming essential tools in semiconductor engineering, particularly for accurately estimating stress distributions.These models are extremely useful tools throughout the design phase, making it easier to create semiconductor devices that are more durable and dependable.Engineers may optimize manufacturing processes by integrating real-time temperature fluctuations and their resulting mechanical reactions into these models, allowing them to adjust process settings and minimize the likelihood of mechanical failures.This rigorous optimization approach not only guarantees the dependability and efficiency of the ultimate semiconductor products but also makes a substantial contribution to enhancing overall production effectiveness and output.As a result, the inclusion of thermoelastic models is critical for the advancement of semiconductor technology and its various applications.
Lord and Shulman 4 introduced their concept of generalized thermoelasticity with a single relaxation time as an alternative to Boit's coupled thermoelasticity theory 5 .Sherief and Saleh 6,7 developed the theory of thermoelastic diffusion.According to this theoretical framework, both mechanical and thermal waves travel at finite speeds, in contrast to the unlimited speeds observed in the coupled theory.References [8][9][10][11][12][13] provide applications of the generalized thermoelastic theory, which is associated with a single relaxation time.Abd El-Latief et al. 14 recently published two novel methods that use four random walks to solve the 1D thermal shock problem in a thermoelastic half space using surrogate-based global optimization.
Thermoelasticity can significantly affect the electrical properties of a material.The thermoelastic effect can change a semiconductor's carrier concentration based on temperature.This phenomenon can lead to changes in the material's electrical conductivity as well as adjustments in its optical features.Gordon et al. 15 were the first to discover the photothermal approach, where they showed that a laser-based device could cause photothermal blooming, leading to the formation of the photothermal lens.Kreuzer 16 demonstrated the efficacy of photoacoustic spectroscopy for precise analysis using laser light sources.There are two types of photothermal processes that happen in semiconductor materials.The thermoelastic (TE) mechanism uses thermal waves to create elastic vibrations, while the electronic deformation (ED) mechanism uses photo-excited free carriers to make periodic elastic deformation 17 .Comprehensive theoretical analyses 18 have described the mechanisms through coupled plasma, thermal, and elastic wave equations.
Many scholars have studied the uncoupled system consisting of plasma, thermal, and elastic equations.Mandelis 19 conducted a study that examined the phenomenon of coupled thermal and plasma waves, also known as "thermoelectronic" waves.McDonald and Wetsel 20 investigated the thermal and acoustic properties of a nonelectronic material.Stearns and Kino 21 , and Todorovic 22 have examined the system of partially coupled plasma, thermal, and elastic wave equations.Todorovic analyzes the reasons for disregarding the relationship between these equations.Furthermore, Todorovic employs the modified Fourier law to deduce a hyperbolic equation that characterizes the process of heat conduction 23 .The equation accurately predicts the limited speed at which heat waves travel through semiconductors.Researchers have studied the interaction between plasma and elastic waves on a semiconductor plate, also known as plasmaelastic (PE) waves 24 .Additional research built upon the existing groundwork, including the work of Song et al. 25 , who examined problems related to the reflection of light in semiconductors using generalized thermoelastic and plasma models.Othman et al. 26 utilized normal mode analysis to study photothermal waves, whereas Abbas and Aly 27 found analytical solutions for waves generated by concentrated laser beams.Later studies by Abbas 28 and Hobiny and Abbas 29 looked at how plasma, thermal, and elastic waves interact with cylinder-shaped cavities in semiconductor materials.Youssef et al. 30 have introduced the hyperbolic two-temperature model in generalized thermoelasticity.Laser pulses cause semiconductor thermoelastic waves, as studied by Taye et al. 31 .Abo-Dahab and Lotfy 32 investigated the problem of two-temperature

Derivation of the fundamental equations
Let V stands for any volume of semiconductor material that is enclosed by a surface A .Consequently, the rule of conservation of energy for V can be expressed as 6 where N = n − n 0 , and M is the energy added to the system when adding one particle without adding heat or work 39 .The presence of superimposed dots signifies time derivatives, and the summation convention rule is applied as usual.
The equation of motion and the symmetry of the stress tensor are given by Apply Gauss' divergence theorem to Eq. (1), and using Eq. ( 2), we get the pointwise form The entropy equation is given by 40,41 Equation (3) after using Eq. ( 4), becomes We now define the Helmholtz free energy function ψ as Eliminating U between Eq. ( 5) and Eq. ( 6), we get The function ψ (as well as all other considered functions) may be defined in terms of the independent vari- ables e ij , T, and N, then by using the chain rule, we obtain Comparing both Eqs.(7) and (8), yields (1) The function ψ is now expanded as a power series of the mentioned independent variables in the given form 42 where θ = T − T 0 , and θ T 0 ≪ 1 .Both a and b are physical constants with dimensions [ J/K ] and [ J.m 3 ].Now, in the medium's natural state, we have Hence, Therefore, we can rewrite Eq. ( 12) as follows (ignoring terms of higher order): Using Eqs. ( 9)- (10) into Eq.( 15), we obtain For the isotropic case, we have the following relations.
where α n represents the discrepancy in deformation potential between the conduction and valence bands 22 .
Furthermore, the parameters a and b can be written as follows: 23 By substituting Eqs. ( 19)- (20) into Eqs.( 16)- (17), we obtain their form in the isotropic case.
By substituting Eq. ( 23) into Eq.( 2) with the definition of the strain tensor e ij = 1 2 (u i,j + u j,i ) , we get We can rewrite Eq. ( 4) in a linearized form as follows: Using the above equation into Eq.( 25), we obtain We now assume the generalized Fourier law of heat conduction as www.nature.com/scientificreports/Next, we apply the divergence operator to both sides of Eq. ( 29) and use Eq. ( 28) along with its time derivative to arrive at the equation of heat conduction.
For low electric field, the current carrier density vector J is given by 43 where the bold letters refer to vector quantities, and the magnetic induction B is given by the relation Equation ( 31) can be rewritten for a semi-conductor medium in terms of chemical potential as follows 44 We postulate an analogous equation for the current density vector J , similar to Eq. ( 29) for the heat flux vector The significance of the aforementioned modification will be addressed in section "Particle diffusion relaxation time : significance and impact on current carrier" at a later point.
The conservation of charge equation for semiconductors is 45 By applying the divergence operator to both sides of Eq. ( 34) and utilizing Eq. ( 24) and Eq. ( 35), we obtain Therefore, our model is formulated, comprising three equations: the equation of motion (Eq.( 26)), the equation of heat conduction (Eq.( 30)), and the particle diffusion equation (Eq.( 36)), which are supplemented by the constitutive equations (Eqs.( 23)-( 24)).The body force term in the equation of motion can be attributed to the Lorentz force, which is expressed as.
Lorentz's force explains the mathematical equations and the physical significance of the forces exerted on charged particles when they traverse through a region of space that contains both an electric field and a magnetic field.

Formulation of the problem
Consider a homogeneous isotropic semiconductor material occupying the region of the half space 0 ≤ x < ∞ that is taken to be traction free and is subjected to a time dependent thermal shock where the temperature is assumed to be a known function of time on the surface of the half plane.So, due to the physics of the problem, all the functions will depend on x and t only; therefore, the displacement components become u x = u(x, t) , u y = u z = 0.
In the absence of body force, heat source, carrier source, and external applied magnetic field.The governing equations can be expressed as In addition to the constitutive equations ( 30) (41) T h e n o n -d i m e n s i o n a l v a r i a b l e i s g i v e n b y t h e f o l l o w i n g r e l a t i o n s 38)- (39) in the dimensionless form become where The initial conditions of the problem are taken to be homogeneous w hile the boundary conditions are assumed to be where H(t) is the Heaviside unit step function defined by

Solution in the Laplace transform domain
We now define the Laplace transform for a function f (t) by the relation Applying Laplace transform on Eqs. ( 44)- (45), and using the homogeneous initial conditions, we get where D = d dx .By simultaneously solving Eqs. ( 54)-( 55), we obtain the following differential equation.
where The bounded solutions of Eq. ( 59) as x → ∞ are given by where A i (s), B i (s) and D i (s) are parameters to be determined and k 1 , k 2 and k 3 are the roots with positive real part of the following characteristic equation The roots k 1 , k 2 , and k 3 of the above characteristic equation are given by where p = ξ 2 1 − 3ξ 2 , q = 1 3 sin −1 (r) and r = − . By substituting Eqs. ( 60)-(61) into Eqs.( 44)- (45), after some manipulations, we obtain (59)

Numerical result and discussion
In order to present numerical results in an effective way, a Silicon material was chosen for purposes of numerical evaluation, whose parameters in SI units are listed in Table 1 [46][47][48][49][50] .
In order to invert the Laplace transform in the above Eqs.( 64)-( 65), we adopt a numerical inversion method based on a Fourier series expansion 51 .The numerical inversion method used to find the solution in the physical domain is listed in [52][53][54][55][56][57][58] .The FORTRAN programming language generated the numerical code.The MATLAB software, with its powerful graphical capabilities, visually represented the numerical results, maintaining a precision of ten digits in the numerical algorithm.
Figures 1, 2, 3 show that temperature, displacement, and stress act in line with the generalized theory of thermoelasticity, which predicts a finite speed for small times, and behave like the coupled theory, which predicts an infinite speed for large times.
Figure 1 is plotted to show the variation of dimensionless temperature θ(x, t) against x−axis for three different instants of time, namely t = 0.1, 0.2 and t = 0.3 .We observe that at any fixed time, the temperature distribution reaches its maximum value at the boundary, aligning with the thermal boundary condition that subjects the half space's boundary surface to a thermal shock.Inside the medium, the values slowly decay.For small values of time (t = 0.10, say) , the values drop to zero, and the sharp temperature decreases to zero due to the effect of the thermal wave.This behavior is consistent with the generalized theory of thermoelasticity, which predicts finite wave speeds for short durations.Over longer values of time ( t = 0.2or0.3) ,the temperature profile aligns with the coupled theory, exhibiting characteristics of infinite speed propagation.This shift indicates that initially, thermal effects are localized, but over longer times, the effects spread rapidly throughout the material.
Figure 2 depicts the variation of dimensionless displacement u(x, t) against x− axis for three different instants of time, namely t = 0.1, 0.2 and t = 0.3 .We notice that for any fixed value of time, the displacement records a  negative value at the boundary.This initial value indicates a significant mechanical response to the thermal input.
The displacement value then increases with distance, reaching its maximum positive value, before gradually decreasing to zero.In addition, when the values of time increase, the magnitude of the displacement increases on the whole domain and cut x-axis more lately, signifying an increase in the rate of mechanical deformation.We plot Fig. 3 to illustrate the variation of the dimensionless stress component σ xx (x, t) with distance x at various value of time.Initially, the stress distribution shows a sharp gradient near the surface, indicating significant stress generation due to thermal effects.For any fixed value of time, the magnitude of stress increases from zero values at the boundary to the maximum values, then decreases with the increase in distance, and finally diminishes to zero.This two-way interaction between the thermal and mechanical waves in the non-dimensional temperature distribution and non-dimensional stress distribution shows how complicated it is for silicon's thermal and mechanical responses to interact with each other.Also, for small values of time, the locations of the thermal and mechanical wave fronts coincide in Figs. 1 and 3.This is evidence of the accuracy of the numerical calculations provided.Figure 4 illustrates the non-dimensional distribution of the carrier density N(x, t) against x-axis for three different instants of time, namely t = 0.1, 0.2 and t = 0.3 .At the boundary x = 0 , there is a high concentra- tion of particles near the surface, which decreases as they diffuse deeper into the silicon medium.Raising the temperature leads to an augmentation in the movement of electrons from the valence band to the conduction band, resulting in an improvement in conductivity and a decrease in resistance.This behavior is in line with the principles of semiconductor physics, where an increase in temperature has a substantial impact on the number of free charge carriers.This phenomenon demonstrates the impact of temperature on the movement of carriers.
Figure 5 shows the non-dimensional distribution of current carriers I(x, t) against x−axis for three different instants of time, namely t = 0.1, 0.2 and t = 0.3 .At the boundary, there is a large value of current carriers on the surface, which decreases as they move into the medium.The initial value represents the prompt reaction of the current carriers to thermal excitation.Over time, the carriers infiltrate more into the silicon, causing a decrease in the original value, which suggests that the carriers are spreading throughout the material.
Figure 6 illustrates the distribution of non-dimensional electrochemical potential energy M(x, t) for different values of time.The surface of the material has the largest potential energy, which gradually diminishes to zero  Based on the Fermi-Dirac distribution, the chemical potential-which stands for the energy required to get an electron into the system-determines the probability that energy levels will be occupied.The chemical potential has a direct impact on the carrier concentration, which consists of holes in the valence band and electrons in the conduction band.To check the validation of the results obtained in Figs. 5, and 6, we can represent the relation between chemical potential and the current carriers for different values of time (see, Fig. 7). Figure 7 depicts the correlation between non-dimensional current carriers I(x, t) and electro-chemical potential energy M(x, t) at different values of time in semiconductor materials.The distribution illustrates the interdependence between the current carriers and electrochemical potential, where changes in one entity impact the other 59 .
Based on Figs. 8, 9, 10, it is seen that the profile of the exhibited functions begins after a period of delay, during which the wave arrived to this deep inside the medium.As the value of x increases, the wave starts later and propagates farther into the media.This phenomenon is consistent with the physical reality, as the waves penetrate farther into the medium and the wave front eventually reaches a greater depth as the elapsed time grows.Furthermore, the function's value has a greater magnitude for shorter distances (smaller x ) as opposed to deeper distances.Figure 8 illustrates the non-dimensional current density N(x, t) at various depths within the silicon medium.The concentration of particles decreases as depth increases, indicating the existence of the diffusion process.The phenomenon of wave arrival delay at increasing depths is in accordance with the physical phenomenon of wave propagation.Analogous to the dispersal of particles, the concentration of current carriers diminishes as depth and duration increase, as shown in Fig. 9.By observing Figs. 8, 9, it is evident that the high concentration on the boundary surface diminishes as carriers spread out into the medium by diffusion.This distribution supports the theoretical model of charge carrier dynamics in semiconductors, where thermal stimulation causes a progressive movement of carriers into deeper regions.
Figure 10 demonstrates that the distribution of electrochemical potential energy grows with time at all depths, suggesting a progressive accumulation of electrochemical potential energy inside the medium.As the depth within the medium increases, the electro-chemical potential energy decreases.This behavior can be attributed to diffusion and transport mechanisms, in which the distribution is driven by a gradient of potential energy and eventually reaches a stable state.

Consent to participate
All authors consent to participate to this publication.

Particle diffusion relaxation time τ 1 : significance and impact on current carrier
In ideal circumstances, the current carriers flowing through semiconductors should have a continuous function.This assumption simplifies theoretical models and circuit analysis, allowing us to apply well-established principles like Ohm's law and Kirchoff 's laws.However, under certain practical scenarios, the current carriers could appear discontinuous or contain rapid transitions.Such situations arise from inherent limitations in modeling assumptions, manufacturing imperfections, or complex interactions occurring in advanced electronic systems.So, it was necessary to address the limitations of theoretical models in explaining this phenomenon in our model.Consequently, the modification of Ohm's law, as described in Eq. (34), is needed.To study the effect of this modification according to different values of particle diffusion relaxation time, a Video S1 displays the current carrier distribution I(x, t) inside the medium, which τ 1 varies from 0 to 0.02 .We note that for τ 1 = 0, the current carriers behave as a smooth continuous function consistent with the theoretical models.However, for τ 1 > 0 , the current carriers suffer a sudden jump in their value at the wavefront location, which becomes more noticeable as τ 1 increases.This covers several experimental situations, where τ 1 is an arbitrary constant that takes special values according to different physical circumstances.

Validation of our mathematical model in comparison with many earlier models.
In this section, we will present a generalized mathematical model based on specific mathematical parameters ζ i , i = 1,2, . . ., 8 .We can use this model to derive every earlier model 22,23,60 as a special case, which are included in Table 2.The proposed governing equations are as follows: where γ u is the elastic coupling factor.By substituting the non-dimensional variables from the previous Eq.( 43), we obtain where ω 1 = E g β 1 τ R C E β 2 η( +2µ) , and α 8 = γ u q 2 0 η .We obtained the solution for the various previous models using the same methods in section "Formulation of the problem", but we don't present them here to avoid expanding the article's body.The Tables 3, 4, 5 are critical for verifying our model's precision.They allow us to compare our model (I) with other models (II, III) by examining non-dimensional temperature, stress and current carrier distributions at various locations ( x ) and two instants of time (t = 0.05, t = 0.10).These comparisons emphasize the model's ability to accurately forecast temperature and stress behavior.
(69) The presented model  11a, b), nondimensional stress distribution (see Fig. 11c, d), nondimensional displacement (see Fig. 12a, b), and nondimensional number of particles (as in Fig. 12c, d).In those figures, labels (a) and (c) indicate the case of small time ( t = 0.05 ), while those labeled by (b) and (d) indicate the case of large time ( t = 1 ).The non-dimension relaxation parameters for the generalized model III are taken to be τ 0 = τ 1 = 0.02 , while those values are τ 0 = 0.02, τ 1 = 0.01 for our model I.
From these figures, we note that (i) The temperature distribution for small values of time in our model has a finite speed, while the others predict an infinite speed, and the wave front location appears at x = 0.35 .This advantage distinguishes our model, while the other two models behave like the coupled theory of thermoelasticity.Therefore, our model eliminates the paradox in the other two models.The temperature distributions for all models are identical for large values of time.(ii) The normal stresses component has the same mechanical wave front location for all models at any time, while the thermal wave front location appears only in our model for a short time.There is a slight difference between all models in the tensile region, while it is more prominent in the compressive (iii) The displacement component u distribution has a small change between all models in the case of small time, while there is a prominent change in the case of large time.(iv) The number of particle distribution N in our model records a smaller value rather than other models and has rapidly decay to zero.(v) All functions field in our model have a smaller value than the other models and rapid decay in the semiconductor media.

Conclusion
This study presents an innovative mathematical framework for examining the characteristics of semiconductor elastic materials when subjected to an external magnetic field.The main goal is to provide a comprehensive mathematical model that encompasses the interplay of plasma, thermal, and elastic waves in semiconductor materials.The mathematical modeling of semiconductors now includes the integration of electrochemical potential, marking a significant development.The electrochemical potential in semiconductors is vital as it controls the behavior and concentration of charge carriers, hence directly impacting the performance and functioning

Figure 1 .
Figure 1.Non dimensional temperature distribution θ(x, t) at different values of time.

Figure 2 .
Figure 2. Non dimensional displacement distribution u(x, t) at different values of time.

Figure 3 .
Figure 3. Non dimensional normal stress distribution σ xx (x, t) at different values of time.

Figure 4 .
Figure 4. Non dimensional carrier density distribution N(x, t) at different values of time.

Figure 5 .
Figure 5. Non dimensional current carrier distribution I(x, t) for different values of time.

Figure 6 .Figure 7 .
Figure 6.Non dimensional electro-chemical potential energy distribution M(x, t) for different values of time.

Figure 8 .
Figure 8. Non dimensional number of particle distribution N(x, t) for different depths inside the medium.

Figure 9 .
Figure 9. Non dimensional current carrier distribution I(x, t) for different depths inside the medium.

Figure 10 .
Figure 10.Non dimensional electro-chemical potential energy distribution M(x, t) for different depths inside the medium.

Figure 11 .
Figure 11.Comparison between the temperature and Normal stress distributions for different models at small and large instant of times.

Figure 12 .
Figure 12.Comparison between the displacement and number of particles distributions for different models at small and large instant of times.

Table 1 .
Physical material parameters of silicon in SI units at T 0 = 293K.

Table 2 .
Limiting cases of the suggested model with previous models.

Table 3 .
Non-dimensional temperature distributions at various locations ( x ) and two instants of time (t = 0.05, t = 0.10).Figures 11, 12 represent a comparison between the above different models (Model I, black line), (Model II, blue line), and (Model III, red line) for non-dimensional physical variable fields, namely, nondimensional temperature distribution (see Fig.

Table 4 .
Non-dimensional temperature distributions at various locations ( x ) and two instants of time (t = 0.05, t = 0.10).

Table 5 .
Non-dimensional temperature distributions at various locations ( x ) and two instants of time (t = 0.05, t = 0.10).